# 自适应积分
# Program 5.2 Adaptive Quadrature
# Computes approximation to definite integral
# Inputs: Matlab function f, interval [a0,b0],
# error tolerance tol0
# Output: approximate definite integral

import numpy as np
from numpy import ndarray

np.set_printoptions(formatter={'float': '{:<#6.5g}'.format})


def trap(f, a, b):
    s = (f(a) + f(b)) * (b - a) / 2.0
    return s


def adapquad(f, a0, b0, tol0):
    # allocate memory
    arrL: int = 3000
    a: ndarray = np.empty(arrL, dtype=float)
    b: ndarray = np.empty(arrL, dtype=float)
    tol: ndarray = np.empty(arrL, dtype=float)
    app: ndarray = np.empty(arrL, dtype=float)
    # initialize
    integral: float = 0
    n: int = 1
    stateL: int = 0  # 当前区间划分长度
    calcL: int = 0  # 计算次数
    subL: int = 0  # 划分子区间长度
    a[0] = a0  # 积分小区间 起点 list
    b[0] = b0  # 积分小区间 终点 list
    tol[0] = tol0  # 每个小区间的 积分容差,
    app[0] = trap(f, a0, b0)
    while n > 0:  # n is current 积分区间长度
        iB: int = n - 1  # 最大索引
        c: float = (a[iB] + b[iB]) / 2.0
        old_app = app[iB]
        app[iB] = trap(f, a[iB], c)
        app[iB + 1] = trap(f, c, b[iB])
        calcL += 2
        if np.abs(old_app - (app[iB] + app[iB + 1])) < 3 * tol[iB]:
            integral += app[iB] + app[iB + 1]  # success
            n -= 1  # done with interval
        else:  # divide into two intervals
            b[iB + 1] = b[iB]
            b[iB] = c
            a[iB + 1] = c
            tol[iB] /= 2.0
            tol[iB + 1] = tol[iB]
            n += 1  # go to end of list, repeat
            subL += 2  # 增加两个新的 子区间
        stateL = np.amax((n, stateL))
    #-------- debug
    print(f'state length        : {stateL}')
    print(f'call f(x) times     : {calcL}')
    print(f'sub interval num    : {subL}')
    # ---------------
    a.resize(stateL)
    b.resize(stateL)
    tol.resize(stateL)
    app.resize(stateL)
    integral.resize(stateL)
    print(f'a is\n{a}\n')
    print(f'b is\n{b}\n')
    print(f'tol is\n{tol}\n')
    print(f'app is\n{app}\n')
    print(f'integral is\n{integral}\n')

    return integral


def test():
    f = lambda x: 1 + np.sin(np.exp(3.0 * x))
    integral = adapquad(f, -1, 1, 0.005)
    print(f'adaptive quad is:\n{integral}')


if __name__ == '__main__':
    test()
